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Disruptions in the normal rhythmic functioning of the heart, termed as arrhythmia, often result 
from qualitative changes in the excitation dynamics of the organ. The transitions between different 
types of arrhythmia are accompanied by alterations in the spatiotemporal pattern of electrical 
activity that can be measured by observing the time-intervals between successive excitations of 
different regions of the cardiac tissue. Using biophysically detailed models of cardiac activity we 
show that the distribution of these time-intervals exhibit a systematic change in their skewness 
during such dynamical transitions. Further, the leading digits of the normalized intervals appear 
to 6t Benford’s law better at these transition points. This raises the possibility of using these 
observations to design a clinical indicator for identifying changes in the nature of arrhythmia. More 
importantly, our results reveal an intriguing relation between the changing skewness of a distribution 
and its agreement with Benford’s law, both of which have been independently proposed earlier as 
indicators of regime shift in dynamical systems. 

PACS numbers: 87.19.Hh,05.45.-a 

I. INTRODUCTION 

Many vital physiological processes are characterized by 
rhythmic activity, ranging from the circadian clock reg¬ 
ulating the daily sleep-wake cycle to temporal patterns 
of respiration that occur over a scale of seconds ll|. The 
periodic beating of the heart, that results in constant cir¬ 
culation of oxygenated blood throughout the body, is one 
of the most important of such naturally occurring oscilla¬ 
tory phenomena in the body 0 . Certain types of distur¬ 
bances in the cardiac rhythmicity, referred to as arrhyth¬ 
mia, can severely impair the normal functioning of the 
heart and in the most critical instances, can result in sud¬ 
den cardiac death Q. Such “dynamical diseases” ii, 
i.e., diseases resulting from abnormal activity in an other¬ 
wise intact physiological system, are a significant public 
health burden in developed countries. For example, in 
the United States, diseases of the heart constitute the 
leading cause of death (responsible for about 25% of all 
deaths), of which more than half can be classified as 
sudden cardiac deaths [^-Q. Even in developing coun¬ 
tries, in recent times heart disease has overtaken other 
causes of death, e.g., sudden cardiac deaths contributed 
to about 10% of overall mortality in certain regions in In¬ 
dia, accounting for upto half of all cardiovascular-related 
deaths P. [l^. 

Several studies have shown that early detection of on¬ 
set of arrhythmia resulting in prompt therapeutic inter¬ 
vention significantly improves the chances of surviving 
such episodes [ll|, [l^- Thus, developing methods for 
identifying signs of impending arrhythmic events with 
potentially serious consequences can significantly con¬ 
tribute towards reducing the mortality rate due to sud¬ 
den cardiac death. With this aim in view there have 


been a number of attempts at applying time-series anal¬ 
ysis methods on cardiac activity data in order to extract 
robust indicators of imminent instances of temporal irreg¬ 
ularities in the heart. However, the complexity of heart 
rate dynamics makes it difficult to characterize and dis¬ 
tinguish the temporal signatures of a healthy heart from 
a diseased one (la - fl^ . 

Most studies of cardiac time-series have focused on 
heart rate variability as measured by temporal changes 
in the R-R interval, the duration between successive 
episodes of ventricular depolarization which triggers con¬ 
traction of the lower chambers of the heart. Following the 
pioneering observations connecting decreased variance in 
R-R intervals with higher mortality risk in patients suf¬ 
fering myocardial infarction [lO, El , it is now generally 
accepted that healthy individuals have higher heart-rate 
variability compared to those with diseased hearts (^ . 
However, certain pathological conditions including car¬ 
diac arrhythmia are seen to be extremely irregular [^ . 
In fact, it has been observed that a transition from tachy¬ 
cardia, i.e., abnormally rapid heart-rate, to fibrillation, 
characterized by erratic muscle activity that prevents the 
heart from pumping blood, is marked by a switch from 
relatively more periodic activity to a highly irregular dy¬ 
namical state [^. While the R-R intervals in normal 
sinus rhythm appear to have almost as unpredictable a 
nature as that seen during fibrillation (^ . it has been 
suggested that the “chaoticity” during normal cardiac ac¬ 
tivity arises through interaction of the heart with the ner¬ 
vous system [2^. In contrast, the spatiotemporal chaos 
associated with fibrillation arises from intrinsic instabil¬ 
ities in cardiac excitation dynamics (E| . 

One possible route from tachycardia to hbrillation that 
has been established through extensive simulation studies 
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of models of cardiac electrical activity is the degeneration 
of reentrant spiral wave (corresponding to tachycardia) to 
disordered, turbulent activit y (c haracterizing fibrillation) 
through spiral breakup |28l . l29l| . This dynamical transi¬ 
tion has been reproduced in a wide range of systems, from 
simple, excitable media to biologically detailed models, 
underlining the robustness of the scenario [s^. Thus, the 
study of spatiotemporal dynamics in models of electrical 
activity in cardiac tissue provides another perspective to 
identify indicators for an impending onset of possibly life- 
threatening arrhythmia. 

In this paper, we focus on analyzing time-series data 
obtained from spatially extended models of cardiac ven¬ 
tricular activity in which, by tuning specific physiological 
parameters, one can observe transitions to different dy¬ 
namical regimes representing various classes of arrhyth¬ 
mia. This allows us to look for statistical signatures 
that can help in early detection of arrhythmic episodes, 
where the observed patterns are exclusively due to abnor¬ 
mal excitation activity that characterizes such arrhyth¬ 
mia and unrelated to heart rate variability that arises 
from the influence of the nervous system on the sinus 
node, the natural pacemaker of the heart [Slj. This 
study, therefore, provides a benchmark against which 
analysis of ECG data obtained from clinical studies can 
be compared, enabling distinction of statistical features 
of arrhythmic time-series that are intrinsic to the dy¬ 
namics of heart muscle from those that are a result of 
changes in the autonomic modulation of cardiovascular 
function (achieved through dynamical balance between 
sympathetic and parasympathetic effects (lH). As sig¬ 
nature patterns (if they exist) that indicate transitions 
from one dynamical regime of cardiac activity to another 
may be masked by other effects in reality, establishing 
them through analysis of the model output will allow us 
to look for them in data obtained from experimental or 
clinical studies. 

Here, for our statistical analysis, we have focused on 
the sequence of time-intervals T between successive ex¬ 
citations of ventricular muscle cells (analogous to the 
R-R interval for ECGs) as a representative feature of 
heart rate dynamics (Eig. [T|). An important result of 
our study is that the distribution of these intervals ex¬ 
hibit clearly observable changes in their moments - in 
particular, the skewness - around the onset of qualita¬ 
tively distinct dynamical behavior characterizing various 
arrhythmic episodes (represented by the different pan¬ 
els in Fig. [T|). Intriguingly, we also observe that at these 
transitions, the distribution of the time-intervals appears 
to agree better with Benford’s law (BL), an empirically 
established feature of the frequency distribution of lead¬ 
ing digits of numbers occurring in many phenomena in 
various physical, biological and social contexts [ssl - l^ . 
Both variation in skewness [13, and closer agreement 
with BL have independently been suggested as in¬ 
dicators of regime shifts or phase transitions in different 
systems. Our work not only finds both of these signa¬ 
tures to be indicative of the onset of arrhythmic behav- 
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FIG. 1: Time-series of the transmembrane potential V rep¬ 
resenting local excitation activity in a two-dimensional LRl 
model {L = 400) for different values of the maximum Ca^^ 
channel conductance Gsi, viz., (A) 0.005, (B) 0.04 and (C) 
0.065 mS cm~^. The distinct nature of the corresponding spa¬ 
tiotemporal dynamics, viz., rigid rotation of a spiral similar 
to monomorphic tachycardia (A), chaotic meandering of spi¬ 
ral core representing polymorphic tachycardia (B), and spa¬ 
tiotemporal chaos indicating fibrillation (C), is visually ap¬ 
parent in the varying pattern of intervals between successive 
excitations. It is quantified in terms of the sequence of time- 
intervals {Tn} (n = 1,2,...) between each pair of consecutive 
local supra-threshold depolarizations (two such intervals are 
shown in panel C using double-headed arrows). We consider a 
local supra-threshold excitation event to have occurred when 
V exceeds —50 mV (broken line). 


ior, but additionally suggests that these two phenomena 
(viz., increasing skewness and agreement with BL) may 
be related. 


II. MATERIALS AND METHODS 


Model. To simulate spatiotemporal excitation activ¬ 
ity in cardiac muscle under different physiological condi¬ 
tions, we have used a two-dimensional model of ventric¬ 
ular tissue having the generic form 




lionjV, di) ^ 
Cm 


( 1 ) 


where V (mV) is the potential difference across a cellu¬ 
lar membrane, Cm{= l/rFcm“^) is the transmembrane 
capacitance, D is the diffusion constant ( = O.OOlcm^s"^ 
for the results reported in the paper), /io„(/j,Acm“^) is 
the total current density through ion channels on the cel¬ 
lular membrane, and Qi describes the dynamics of gating 
variables of different ion channels. The specific functional 
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form for lion that we have focused on here is that of the 
Luo-Rudy I (LRl) model which describes the ionic cur¬ 
rents in a guinea pig ventricular cell (4?)j |: 

lion — ^Na ^Kl “t“ Ikp ^si “t“ df,, 

where Ino = GNa'm-^hj{V — 54.4) is the fast inward 
Na+ current, I si = Gsidf{V — Esi) is the slow inward 
Ca^+ current where Egi = 7.7 — 13.0287 In ([Ca^+]i) is 
the reversal potential, dependent on the intracellular ion 
concentration [Ca^+], Ik = Gkxxi{V + 77.62), Iki = 
GKiRrioo(t^ + 87.95) and J^p = 0.0183R'p(R-b87.95) are 
three different types of K'*' current, and h = 0.03921 (R-|- 
59.87) is a background current. The currents are deter¬ 
mined by ion channel gating variables m, h, j, d, f and 
X, whose time evolution is described by ordinary differ¬ 
ential equations of the form, de/dt = (coo — Gj/r^, where 
Coo is the steady state value of e (=m, h, d, / and x) 
and Te is the corresponding time constant obtained by 
fitting experimental data. Parameter values used are as 
in Ref. [^, except for Gk which is set to 0.705 mS//rF 
and Gsi that is varied to alter the stability of spiral wave 
dynamics [4l|. We have explicitly verified that our re¬ 
sults are not sensitively dependent on model-specific de¬ 
tails, viz., the description of ion-channel dynamics, by ob¬ 
serving qualitatively similar behavior with another func¬ 
tional form for lion described in the ten Tusscher-Panfilov 
(TP06) model of a human ventricular cell (i^ . 

For numerical simulations, the two-dimensional spa¬ 
tially extended system is discretized on a grid of size LxL 
(= 400 for LRl model and 1024 for TP06 model) with a 
space step of Sx (= 0.0225 cm for LRl model and 0.025 
cm for TP06 model). The equations are solved using a 
forward Euler method with time step St (= 0.05 ms for 
LRl model and 0.02 ms for TP06 model) and a standard 
5-point stencil for the Laplacian describing the spatial 
coupling between the lattice elements. No-flux boundary 
conditions are implemented at the edges. The initial spi¬ 
ral wave state is obtained by generating a broken wave 
front which then dynamically evolves into a curved ro¬ 
tating wave. The movement of the spiral wave core is 
obtained by tracing the trajectory of intersection points 
of iso-contour lines for a pair of dynamical variables of 
the model, viz., V and h in the LRl model [43 |. 

Inter-spike interval time-series. We analyze the 
statistical properties of the time-intervals between suc¬ 
cessive excitations (i.e., depolarization) at specific loca¬ 
tions in the simulated cardiac tissue. Each point in the 
simulation domain is considered to be excited if the corre¬ 
sponding transmembrane potential V crosses a threshold 
value (set equal to —50 mV here, although our results are 
robust with respect to the choice of threshold) from be¬ 
low, i.e., from a hyperpolarized state. The time-interval 
T between two such successive crossings of the thresh¬ 
old is recorded for constructing the data-set, values be¬ 
ing sampled from 400 equally spaced points arranged in 
a regular grid on the simulation domain (for the LRl 
model). Typically time series of 5 s total duration were 
used for our analysis. From the data-sets obtained at dif¬ 


ferent values of Gsi, the corresponding distributions for 
T are obtained and the moments calculated, including 
mean /r, standard deviation a, and skewness, the latter 
being measured by the Pearson’s moment coefficient of 
skewness defined as 7 = E[{X — //)/cr)^]. Using other 
measures for the skewness did not qualitatively alter our 
results. The time-interval distributions obtained for dif¬ 
ferent parameters are also tested for the degree of agree¬ 
ment with Benford’s Law. 

BL and Benford distribution. Named after the 
American physicist F. Benford who made the first-digit 
phenomenon widely known, BL had been noticed in num¬ 
bers associated with a variety of natural (as well as social) 
phenomena as far back as in 1881 by the astronomer S. 
Newcomb. According to this empirical law, numbers be¬ 
ginning with 1 or 2 occur more often than those beginning 
with 8 or 9 [1^. Specifically, the probability of the first 
or leading digit of such numbers being i (i = 1, 2,... 9) 
is given by the Benford distribution: 

P{i) = logio + 7 ) • 

The reason for the ubiquity of this distribution has been 
connected to its scale-invariance and base-invariance (4^ . 
Thus, if indeed there is a universal principle underlying 
the distribution of the leading digits of numbers which 
is independent of the units in which the numbers are 
measured or the number base used, then the BL fol¬ 
lows. Simple mathematical arguments have been used 
to show that any distribution of numbers arising from 
natural processes that spans several orders of magnitude 
and is reasonably smooth will obey BL [dj. The Ben¬ 
ford distribution, as mentioned earlier, is seen in many 
empirical data-sets, including those arising in a biological 
context, such as, the distribution of open reading frames 
in prokaryotic and eukaryotic genomes [ 3 ^ . Dynamical 
systems, such as those describing the molecular dynam¬ 
ics of fluids or certain chaotic systems, also exhibit BL 
in the numbers exp ressing coordinates of the generated 
trajectories [4^ . More recently, BL has been used as 

a signature for detecting phase transitions in a quantum 
system [s^. 

In order to compare the distribution of the intervals T 
between successive excitations with that expected from 
BL, we first obtain a set of normalized time-intervals 
t by subtracting the minimum value of the series from 
all intervals T and dividing by the range, i.e., t = 
[T—min(T)]/[max(T)—min(T)]. The leading digits i of 
the normalized intervals t are then extracted as i = 
[|t|/ 10 h°Sio(ld)lJ^ where \z\ and [z\ indicate the abso¬ 
lute value of z and the largest integer not greater than 
2 ; respectively. The distribution of i is then tested for 
agreement with BL using statistical tests for goodness of 
fit. 

Statistical tests for goodness of fit with BL. 

The goodness of fit between the two distributions (the 
empirical and that predicted by BL) is measured by a 
two-sample Kolmogorov-Smirnov test. We have used the 
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function kstest2 implemented in MATLAB which returns 
a test decision for the null hypothesis that both the sets 
are from the same distribution, along with a p-value and 
the KS test statistic k describing the degree of deviation 
from BL. In our study the test statistic is the comparison 
parameter 


k = uia.x{Fc{i) - Bc{i)), 

i 

which measures the maximum distance between the two 
cumulative distributions, F^i) and Bc{i), of the leading 
digits i of normalized time intervals and that expected 
from BL, respectively. A lower value of k implies closer 
agreement with the Benford distribution. 

Apart from the KS test, we have also used the Pear¬ 
son’s chi-squared test to confirm the compliance of the 
empirical distributions with BL. The test statistic 
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quantifies the total magnitude (over all n entries of the 
empirical time-series) of the difference between the two 
probability distributions, F and B, for the leading digits 
of the normalized time intervals and that expected from 
BL, respectively. As for the KS test, a lower value of 
implies closer agreement with the Benford distribution. 


III. RESULTS 

To identify the statistical signatures characterizing dy¬ 
namical transitions to different types of arrhythmia, we 
systematically explore the spatiotemporal dynamics of 
the model systems in different parameter regimes. The 
nature of the excitation activity is varied in a controlled 
manner by changing the kinetic properties of an ion chan¬ 
nel, viz., increasing the maximum Ca^“'" channel conduc¬ 
tance Gsi for the LRl model (keeping all other model 
parameters unchanged) which is known to result in a suc¬ 
cession of dynamical transitions [i^ as seen in Fig.[TJ For 
the TP06 model, the maximum conductance GpCa of the 
sarcolemmal pump current IpCa is increased that even¬ 
tually results in spiral breakup leading to spatiotemporal 
chaos [ 4 ^ . 

Representative images of the spatiotemporal dynamics 
(with LRl model ion-channel kinetics) in the different 
regimes are shown in the top row of Fig. with the tra¬ 
jectory of the spiral core (traced in the first three pan¬ 
els using a light color) exhibiting characteristic changes 
in its qualitative nature. Starting from an initial state 
characterized by a rotating spiral wave, for small values 
of the conductance (e.g. Gsi = 0.005) we observe rigid 
rotation with the core moving around an approximately 
circular trajectory (Fig. [5] A), which corresponds to the 
clinical phenomenon of monomorphic tachycardia. This 
gives way to meandering at higher values of Ggi (— 0.025, 
see Fig. ElB), followed by the appearance of chaotic me¬ 
andering around Gsi = 0.04 (Fig. [5] C) and finally the 


breakup of spiral waves leading to spatiotemporal chaos, 
representative of fibrillation, for values of Gsi > 0.055 
(Fig. El D). 

The panels in the middle row of Fig. E] show the prob¬ 
ability distribution of time intervals between successive 
excitations, T, for the Gsi values corresponding to the 
panels in the top row. We observe that the range of these 
intervals become broader at larger Gsi values as the dy¬ 
namics of the spiral wave becomes more complex. A very 
narrow range of intervals is dominant at low Gsi, as ex¬ 
pected for a rigidly rotating spiral wave having a charac¬ 
teristic period of rotation (Fig. E]E). With increased me¬ 
andering of the core, however, the time interval between 
successive excitations of a local region becomes more ir¬ 
regular, which is manifested as a broader distribution of 
T (Fig.ElF). As the spiral core trajectory becomes even 
more complex, covering a larger portion of the simulation 
domain, we see that the distribution not only widens fur¬ 
ther but also develops multiple peaks at the extremities 
(Fig. El G). Finally, following breakup and spatiotempo¬ 
ral chaos, the time between successive excitations become 
essentially random in character with a distribution that 
spans a relatively large range of T (Fig. El H). 

To see how closely the dynamical process follows the 
Benford distribution in the different regimes, in the bot¬ 
tom row of Fig. El we show the probability distributions 
of the leading digits i of the normalized intervals t. It 
is evident that the distribution of i moves closer to the 
form expected for the Benford distribution (indicated by 
a broken curve) at larger values of Gsi- In fact, the em¬ 
pirical distribution shows the best agreement with BL 
in the spatiotemporally chaotic state corresponding to 
Gsi = 0.065 (Fig. ElL), which is consistent with the cor¬ 
responding time interval distribution being exponential 
in nature - as it is known that values distributed expo¬ 
nentially follow BL. We see that that this distribution of 
leading digits i is closest to BL at the transition points 
corresponding to chaotic meandering {Gsi = 0.040) and 
spiral breakup [Gsi = 0.065). 

To understand the nature of the distributions in the 
different dynamical regimes better, we show how the mo¬ 
ments of the distribution for the time intervals T and that 
for the corresponding leading digits i of the normalized 
intervals vary with increasing Gsi (Fig. El). We observe 
that the mean value of the interval between successive 
excitations steadily rises with Gsi as the complexity of 
the spiral core trajectory increases excepting for a small 
dip around Gsi = 0.055 which is the point of transi¬ 
tion to spiral breakup (Fig. El A). The dispersion of the 
T distribution, measured by its standard deviation ax 
(Fig. El B) shows a similar increasing nature with Gsi 
although, around Gsi = 0.04, where a transition occurs 
from quasiperiodic to chaotic meandering of the spiral 
core, there is a small decrease. The skewness 7 of the 
distribution is the most informative of all the moments 
considered here, as it shows large deviations from zero 
only around critical values of Gsi associated with transi¬ 
tions between different dynamical regimes. In particular, 
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FIG. 2: (A-D) Pseudocolor images of spatiotemporal activity (measured in terms of transmembrane potential V) for the 

two-dimensional LRl model (L = 400) showing the different dynamical regimes obtained by increasing the maximum Ca^+ 
channel conductance Gsi (expressed in units of mS cm“^). The successive panels represent a spiral wave undergoing (A) stable 
rotation (Gsi = 0.005), (B) quasiperiodic meandering (Gai = 0.025) and (C) chaotic meandering (Gsi = 0.04). Further increase 
of Gsi results in breakup of the spiral wave leading to (D) spatiotemporal chaos (Gsi = 0.065). The trajectory of the spiral 
core (the tip of the spiral wave, defined to be a phase singularity) for a duration of 500 ms is indicated in all panels except for 
the one corresponding to chaotic activity where there is a large multiplicity of singularities. (E-H) The probability distribution 
of time intervals T between successive supra-threshold activations of a local region corresponding to the dynamical regimes 
indicated in panels (A-D) respectively. (TL) Probability distribution of the leading digits i of the normalized time intervals 
between successive supra-threshold activations of a local region corresponding to the dynamical regimes indicated in panel 
(A-D) respectively. The broken curve indicates the distribution predicted by Benford’s law. Each distribution in panels E-L is 
obtained by averaging over data collected from many spatial positions in the simulation domain (indicated by points in panel 
A) and also over several realizations, with error bars indicating the standard deviation. 


we notice peaks in the skewness at Gsi = 0.025, 0.04 and 
0.055 which correspond to transition to quasiperiodic me¬ 
andering, chaotic meandering and spiral breakup, respec¬ 
tively (Fig. OC). In order to make the relation between 
the different moments and the dynamical transitions even 
more clear, we have also shown the nature of variation of 
a derived quantity, exp( 77 ’)/(TT, as a function of Gsi- It 
can potentially be used as a statistical indicator for the 
onset of certain types of arrhythmia that may be hard 
to detect by observing the skewness alone. We see from 
Fig. [3] (D) that the measure amplifies the signal indicat¬ 
ing a transition close to Gsi = 0.025 where the spiral 
begins to noticeably meander. 

When we observe the corresponding moments for the 


distribution of leading digits i as a function of Gsi^ we 
note that both the moments /ii and (Fig.[3]E-F) have 
very similar nature of variation, viz., both exhibit dips 
around the values of Gsi at which the different dynami¬ 
cal transitions occur. In contrast, the skewness exhibits 
an almost opposite nature, with peaks occurring at the 
transition points (consistent with increasing skewness of 
the T distribution at these values). This indicates that at 
these points the distribution comes close to the form ex¬ 
pected from BL, as the latter is associated with positively 
skewed distributions. The derived quantity exp('ji)/ai 
conserves this pattern, showing increased values at these 
points. We find that the skewness of T and that of i are 
correlated, the linear correlation coefficient between ■jt 
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FIG. 3: Analysis of various moments for the distributions of the time intervals T between successive supra-threshold activations 
of a local region (A-D) and that of the leading digits i of the normalized time intervals (E-H) as a function of the maximum 
Ca^"*" channel conductance Gsi in the two-dimensional LRl model. Variation in (A) the mean fir, (B) standard deviation (jt, 
(C) skewness 7 t measured in terms of the Pearson’s moment coefficient and (D) the derived quantity exp{'yT)I ctt, correspond 
to the distribution of the time intervals T, while the variation shown for (E) the mean pi, (F) standard deviation ai, (G) 
skewness 7 i (again measured in terms of the Pearson’s moment coefficient), and (H) the quantity exp{'yi)/ai, correspond to the 
distribution of the leading digits i. Large changes in both the skewness measures (7 t and 7 i), and to an extent, the dispersion 
measures {ax and ai) correspond to successive dynamical transitions between rigid rotation, quasiperiodic meander and chaotic 
meander of the spiral core, finally giving rise to spatiotemporal chaos resulting from spiral breakup. The measure exp{‘yi)/ai 
combines the information obtained from the variation of standard deviation and skewness, enabling it to indicate some of the 
dynamical transitions more clearly. The linear correlation coefficient between the skewness of T and that of i is = 0.67 
(p-value= 0.001). 


and 7 i being = 0.67 (p-value= 0.001). This indicates 
an inter-dependence between the variations in skewness 
of the time-interval distribution and that of the leading 
digit distribution, that occur at different dynamical tran¬ 
sitions. 

To quantify how closely the system obeys BL in the 
different dynamical regimes, we show the results of dif¬ 
ferent statistical tests for goodness of fit between the em¬ 
pirical and Benford distributions. Fig 2] (A) shows the 
Kolmogorov-Smirnov (KS) test statistic as a function of 
the Ca^+ channel conductance Gsi which clearly indi¬ 
cates that the distribution of leading digits follow BL 
most closely, indicated by dips in the test statistic, at 
the values of Gsi characterizing the different dynamical 
transitions, viz., Gsi = 0.025, 0.04 and 0.055 (p-values for 
the statistic are effectively zero for all Gsi)- Note that 
low values of the KS test statistic (i.e., better agreement 
with BL) are associated with high positive skewness of 
the distribution of leading digits of the normalized inter¬ 
vals. This is underlined by the strong negative correla¬ 
tion between the test statistic k and the skewness ji (see 
Fig.2]B), having a linear correlation coefficient r = —0.96 


(p= 10“^^). As mentioned earlier, this is consistent with 
the fact that BL is associated with distributions having 
high positive skewness. Fig. |4] (C) shows the result of 
another statistical test, viz., Pearson’s Chi-squared test, 
with the lowest values of Pearson’s error - corresponding 
to better agreement with BL - occurring at the same val¬ 
ues of Gsi where the dynamical transitions occur. The 
points at which the empirical distribution best matches 
BL are seen to be consistent for the two tests (the dips 
of the two test statistics occurring at the same values). 

The change in skewness and the agreement with BL 
around the transition to spiral breakup shown here do 
not appear to be model specihc. We have explicitly veri¬ 
fied that qualitatively similar results can be obtained by 
using the TP06 model for a human ventricular cell. The 
top row of Fig. [S] shows images representative of the spa¬ 
tiotemporal dynamics of a two-dimensional medium in 
different regimes as the maximum value of conductance 
GpCa of the sarcolemmal pump current IpCa is increased. 
For small values of GpCa a single spiral wave is seen to 
rotate stably in the medium, until around GpCa = 0.125 
nS pF“^ it undergoes breakup and degenerates into spa- 
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FIG. 4: (A) Deviation of the distribution of leading digits i of the normalized time intervals t from Benford’s law in the 

two-dimensional LRl model measured in terms of the Kolmogorov-Smirnov test statistic k and shown as a function of the 
maximum Ca^^ channel conductance Gsi- (B) There is a strong negative correlation (r = —0.96, p-value = 10“^^) between k 
and the skewness of the leading digits 7 i. (C) The error in using BL for describing the empirical distribution of leading digits 
i of the normalized time intervals t is measured using Pearson’s chi-squared test and shown as a function of Gsi- Agreement 
with BL improves whenever there is a dynamical transition, as seen by dips in k and the test statistic for values of Gsi 
where successive transitions between rigid rotation, quasiperiodic meandering and chaotic meandering of the spiral core and 
spatiotemporal chaos occur. 



FIG. 5: (A-D) Pseudocolor images of spatiotemporal activity (measured in terms of transmembrane potential V) for the two- 
dimensional TP06 model {L — 1024) showing the different dynamical regimes obtained by increasing the maximum conductance 
GpCa of the sarcolemmal pump current (expressed in units of nS pF”*^). The successive panels represent a single spiral wave at 
GpCa = 0.1 (A), the initiation of spiral breakup at GpCa = 0.125 (B) and successive states leading to spatiotemporal chaos at 
GpCa = 0.15 (C) and GpCa = 0.175 (D). (E-H) Probability distribution of the leading digits i of the normalized time intervals 
between successive supra-threshold activations of a local region corresponding to the dynamical regimes indicated in panel 
(A-D) respectively. The broken curve indicates the distribution according to Benfords law. The analysis shown in panels E-H 
is performed for data obtained from a grid of 225 equally spaced points in the two-dimensional medium. 
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FIG. 6: (A-B) Deviation of the distribution of leading dig¬ 
its i of the normalized time intervals t from BL in the 
two-dimensional TP06 model measured in terms of (A) the 
Kolmogorov-Smirnov test statistic k and (B) Pearson’s Chi- 
squared test statistic shown as a function of the max¬ 
imum conductance GpCa of the sarcolemmal pump current 
(expressed in units of nS pF“^). Agreement with BL is high¬ 
est around the value of GpCa = 0.125 where breakup of the 
spiral wave is initiated, as seen by dips in k and the test 
statistic. (C) Skewness 7 i in the distribution of the leading 
digits of the normalized time interval, with the peak occur¬ 
ring at GpCa = 0.125 corresponding to the spiral breakup 
transition point. (D) There is a strong negative correlation 
(r = —0.99, p-value = 10~^) between k and the skewness of 
the leading digits Xi- 


tiotemporal chaos for higher values of GpCa- The panels 
in the bottom row show the distributions of the leading 
digits for the normalized intervals at the corresponding 
values of GpCa indicating that the fit with Benford distri¬ 
bution is closest during the dynamical transitions. This 
is rigorously shown in Fig. [5] where the results of differ¬ 
ent statistical tests for goodness of fit are given. Panels 
(A) and (B) show that both the KS and Pearson’s Chi- 
squared tests point to a better agreement with BL at the 
transition point GpCa = 0.125 nS pF“^. As with the 
results for LRl model reported earlier, the transition is 
also associated with a large deviation in the skewness 7 i 
(Fig. [51 C) and a strong negative correlation between the 
test statistic k and the skewness (Fig. [51 D) with a linear 
correlation coefficient r = —0.99 {p = 10“®). 

To summarize the results, around the parameter values 
where the transitions between dynamical regimes rep¬ 
resentative of different types of cardiac arrhythmia oc¬ 
cur, we observe both higher positive skewness and closer 
agreement with BL (as indicated by statistical tests). We 
note that both increased skewness and better match with 
BL have independently been suggested earlier as possi¬ 
ble signatures for dynamical transitions, although not in 
the context of physiology or clinical applications. Apart 
from the potential utility of this observation for devising 


robust indicators of the onset of life-threatening distur¬ 
bances in the cardiac rhythm, it suggests a deep relation 
between the appearance of BL in natural phenomena and 
the degree of skewness in the distributions of the under¬ 
lying variables. 


IV. DISCUSSION 

Statistical analysis of data (in particular, ECG) that 
is representative of cardiac functionality can provide us 
with effective signatures for the detection of arrhyth¬ 
mia at an early stage. Despite such analyses, certain 
kinds of arrhythmia fail to get detected merely due to 
the complexity involved. Part of the difficulty lies in 
cardiovascular activity being a joint outcome of intrinsic 
spatiotemporal excitation dynamics in heart muscle and 
modulation of the sinus node by the sympathetic and 
parasympathetic nervous system. Here we have studied 
biophysically detailed models of ventricular activity to 
infer signatures of dynamical transitions characterizing 
the onset of different kinds of arrhythmia. This makes it 
possible to disentangle the effects of intrinsic excitation 
dynamics in the heart from the influence of the nervous 
system. In principle, it allows identification of patterns 
that may alert one to impending harmful disruptions in 
the rhythmic activity of the heart, but which could be 
masked in the ECG signal by autonomic modulation ef¬ 
fects. Our results show that as the spatiotemporal ex¬ 
citations in the ventricles become more disordered, lead¬ 
ing to phenomena identified with sustained monomorphic 
and polymorphic tachycardia as well as the onset of hb- 
rillation, these transitions are marked by characteristic 
changes in statistical moments associated with the dis¬ 
tribution of inter-activation intervals. In addition, the 
leading digits of these intervals show a closer agreement 
with BL at the transition points. Our result can poten¬ 
tially be applied in augmenting algorithms used in im¬ 
planted devices (ICDs) for detecting transitions to pos¬ 
sible life-threatening arrhythmia so as to initiate a pro¬ 
gram of treatment [49|. Thus, when continual monitoring 
of heartbeat time-series shows either increased skewness 
or a closer agreement to BL, it may signal a transition 
in the dynamical state of the heart so that suitable pac¬ 
ing therapy can be started. However, for such an appli¬ 
cation, our observations would need to be validated in 
ECG recordings made during tachycardia and onset of 
fibrillation in live animals. Such validation is necessary 
in view of the limitations of the present study, involv¬ 
ing as it does two-dimensional monodomain models of 
homogeneous cardiac tissue. 

Examining the passage from normal cardiac activity 
to different arrhythmic regimes from the perspective of 
phase transitions can provide novel insights, as has been 
pointed out by several earlier studies. For instance, 
power-law behavior, which characterize critical phenom¬ 
ena in physical systems, have been reported in R-R inter¬ 
val fluctuations and are seen to be remarkable predictors 
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of arrhythmic death, with a steeper negative slope of the 
power spectra (in log scale) clearly distinguishing a dis¬ 
eased heart from a healthy one . More recently, it has 
been shown that phase transition-like dynamics is exhib¬ 
ited by healthy human heart rate, indicated by long range 
correlations which is a hallmark of criticality |5l| . In 
contrast, the dynamics of an abnormal heart rate reveals 
significant digression from critical behavior. In addition, 
scale invariance, which is seen in systems close to criti¬ 
cal point, has been shown to be indicative of a healthy 
heart - with its absence being a statistical feature that 
can alert us about pathological conditions [s^. Our re¬ 
sults provide yet another connection between onset of 
arrhythmia and phase transitions by showing that sharp 
changes in the skewness of the distribution of dynami¬ 
cal variables, that has previously been associated with 
dynamical transitions in other systems [13, , can po¬ 

tentially act as a robust indicator of transitions between 
monomorphic tachycardia, polymorphic tachycardia and 
fibrillation in the heart. 

As mentioned earlier, the appearance of BL has also 
been linked to phase transitions in physical systems [1^ . 
As the Benford distribution follows from re quir ement 
of scale invariance of the underlying numbers |45j |. it is 
tempting to connect this with the scale invariance of dis¬ 
tribution of dynamical variables associated with critical 
points at which transitions occur. We observe from our 
results that parameter regimes that give rise to relatively 
broad distributions of the time-intervals (as is the sit¬ 
uation for spatiotemporally chaotic states) show better 
agreement with BL than those associated with highly 
confined distributions (as in the case of a rigidly rotating 


spiral). However, the occurrence of chaos in dynami¬ 
cal systems by itself does not guarantee that BL will be 
obeyed H^. Thus, it appears that the appearance of 
Benford distribution is more closely associated with the 
onset of dynamical transitions rather than the specific 
nature of the dynamical states on either side of the tran¬ 
sition point. We also note that distributions that follow 
BL are, in general, associated with high skewness [s^. 
This suggests that the highly skewed nature of distribu¬ 
tions during dynamical transitions and the observation 
of better agreement with BL at those points may not 
be independent of each other. Thus, our results imply 
that the increased skewness associated with regime shifts 
and the appearance of Benford distribution during phase 
transitions - which have been reported earlier in different 
contexts - are, in fact, related. 
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